Load libraries

In this step, we are loading the libraries which are required for analysis and plotting of the data.

library(Seurat)
## Attaching SeuratObject
library(ggplot2)
library(patchwork)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(clusterProfiler)
## 
## clusterProfiler v4.2.2  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
## 
## If you use clusterProfiler in published research, please cite:
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141
## 
## Attaching package: 'clusterProfiler'
## The following object is masked from 'package:stats':
## 
##     filter
library(enrichplot)
library(clusterProfiler)
library(enrichplot)
# we use ggplot2 to add x axis labels (ex: ridgeplot)
library(stringr)
library(celldex)
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
## 
##     count
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:clusterProfiler':
## 
##     rename
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## The following object is masked from 'package:clusterProfiler':
## 
##     slice
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## Loading required package: GenomeInfoDb
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
## 
## Attaching package: 'SummarizedExperiment'
## The following object is masked from 'package:SeuratObject':
## 
##     Assays
## The following object is masked from 'package:Seurat':
## 
##     Assays
library(ensembldb)
## Loading required package: GenomicFeatures
## Loading required package: AnnotationDbi
## 
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:clusterProfiler':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: AnnotationFilter
## 
## Attaching package: 'ensembldb'
## The following object is masked from 'package:clusterProfiler':
## 
##     filter
## The following object is masked from 'package:dplyr':
## 
##     filter
## The following object is masked from 'package:stats':
## 
##     filter
library(SingleR)
## 
## Attaching package: 'SingleR'
## The following objects are masked from 'package:celldex':
## 
##     BlueprintEncodeData, DatabaseImmuneCellExpressionData,
##     HumanPrimaryCellAtlasData, ImmGenData, MonacoImmuneData,
##     MouseRNAseqData, NovershternHematopoieticData
library(limma)
## 
## Attaching package: 'limma'
## The following object is masked from 'package:BiocGenerics':
## 
##     plotMA
library(DESeq2)
library(org.Mm.eg.db)
## 
library(rhdf5)
library(hdf5r)
## 
## Attaching package: 'hdf5r'
## The following object is masked from 'package:rhdf5':
## 
##     h5version
## The following object is masked from 'package:SummarizedExperiment':
## 
##     values
## The following object is masked from 'package:GenomicRanges':
## 
##     values
## The following object is masked from 'package:S4Vectors':
## 
##     values
library(pheatmap)
library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.10.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##   genomic data. Bioinformatics 2016.
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## ! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
##    in the original pheatmap() are identically supported in the new function. You 
##    can still use the original function by explicitly calling pheatmap::pheatmap().
## 
## Attaching package: 'ComplexHeatmap'
## The following object is masked from 'package:pheatmap':
## 
##     pheatmap
library(CellChat)
## Loading required package: igraph
## 
## Attaching package: 'igraph'
## The following object is masked from 'package:GenomicRanges':
## 
##     union
## The following object is masked from 'package:IRanges':
## 
##     union
## The following object is masked from 'package:S4Vectors':
## 
##     union
## The following objects are masked from 'package:BiocGenerics':
## 
##     normalize, path, union
## The following object is masked from 'package:clusterProfiler':
## 
##     simplify
## The following objects are masked from 'package:dplyr':
## 
##     as_data_frame, groups, union
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
library(patchwork)
options(stringsAsFactors = FALSE)
library(NMF)
## Loading required package: pkgmaker
## Loading required package: registry
## 
## Attaching package: 'pkgmaker'
## The following object is masked from 'package:S4Vectors':
## 
##     new2
## Loading required package: rngtools
## Loading required package: cluster
## NMF - BioConductor layer [OK] | Shared memory capabilities [NO: bigmemory] | Cores 3/4
##   To enable shared memory capabilities, try: install.extras('
## NMF
## ')
## 
## Attaching package: 'NMF'
## The following objects are masked from 'package:igraph':
## 
##     algorithm, compare
## The following object is masked from 'package:S4Vectors':
## 
##     nrun
library(ggalluvial)

load data

# saline
saline<- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/Saline_sample_round_one/outs/"

list.files(saline)
##  [1] "analysis"                      "cloupe.cloupe"                
##  [3] "filtered_feature_bc_matrix"    "filtered_feature_bc_matrix.h5"
##  [5] "metrics_summary.csv"           "molecule_info.h5"             
##  [7] "possorted_genome_bam.bam"      "possorted_genome_bam.bam.bai" 
##  [9] "probe_set.csv"                 "raw_feature_bc_matrix"        
## [11] "raw_feature_bc_matrix.h5"      "spatial"                      
## [13] "spatial_enrichment.csv"        "web_summary.html"
saline <- Load10X_Spatial(
  saline,
  filename = "filtered_feature_bc_matrix.h5",
  assay = "spatial",
  slice = "saline",
  filter.matrix = TRUE,
  to.upper = FALSE,
  image = NULL)
## Warning in sparseMatrix(i = indices[] + 1, p = indptr[], x = as.numeric(x =
## counts[]), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
# ntm
ntm <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/ntm_sample_round_two/outs/"

list.files(ntm)
##  [1] "analysis"                      "cloupe.cloupe"                
##  [3] "filtered_feature_bc_matrix"    "filtered_feature_bc_matrix.h5"
##  [5] "metrics_summary.csv"           "molecule_info.h5"             
##  [7] "ntm_only_web_summary.html"     "possorted_genome_bam.bam"     
##  [9] "possorted_genome_bam.bam.bai"  "probe_set.csv"                
## [11] "raw_feature_bc_matrix"         "raw_feature_bc_matrix.h5"     
## [13] "spatial"                       "spatial_enrichment.csv"
ntm <- Load10X_Spatial(
 ntm,
  filename = "filtered_feature_bc_matrix.h5",
  assay = "spatial",
  slice = "ntm",
  filter.matrix = TRUE,
  to.upper = FALSE,
  image = NULL)
## Warning in sparseMatrix(i = indices[] + 1, p = indptr[], x = as.numeric(x =
## counts[]), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
# bcg

bcg <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/bcg_sample_round_one/outs/"

list.files(bcg)
##  [1] "analysis"                      "cloupe.cloupe"                
##  [3] "filtered_feature_bc_matrix"    "filtered_feature_bc_matrix.h5"
##  [5] "metrics_summary.csv"           "molecule_info.h5"             
##  [7] "possorted_genome_bam.bam"      "possorted_genome_bam.bam.bai" 
##  [9] "probe_set.csv"                 "raw_feature_bc_matrix"        
## [11] "raw_feature_bc_matrix.h5"      "spatial"                      
## [13] "spatial_enrichment.csv"        "web_summary.html"
bcg <- Load10X_Spatial(
  bcg,
  filename = "filtered_feature_bc_matrix.h5",
  assay = "spatial",
  slice = "bcg",
  filter.matrix = TRUE,
  to.upper = FALSE,
  image = NULL)
## Warning in sparseMatrix(i = indices[] + 1, p = indptr[], x = as.numeric(x =
## counts[]), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
# ntm + bcg
ntm_bcg <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/ntm_bcg_sample_round_two/outs/"

list.files(ntm_bcg)
##  [1] "analysis"                      "cloupe.cloupe"                
##  [3] "filtered_feature_bc_matrix"    "filtered_feature_bc_matrix.h5"
##  [5] "metrics_summary.csv"           "molecule_info.h5"             
##  [7] "ntm+bcg_web_summary.html"      "possorted_genome_bam.bam"     
##  [9] "possorted_genome_bam.bam.bai"  "probe_set.csv"                
## [11] "raw_feature_bc_matrix"         "raw_feature_bc_matrix.h5"     
## [13] "spatial"                       "spatial_enrichment.csv"
ntm_bcg <- Load10X_Spatial(
ntm_bcg,
  filename = "filtered_feature_bc_matrix.h5",
  assay = "spatial",
  slice = "ntmbcg",
  filter.matrix = TRUE,
  to.upper = FALSE,
  image = NULL)
## Warning in sparseMatrix(i = indices[] + 1, p = indptr[], x = as.numeric(x =
## counts[]), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

Data normalization

saline <- SCTransform(saline, assay = "spatial", verbose = FALSE)

ntm <- SCTransform(ntm, assay = "spatial", verbose = FALSE)

bcg <- SCTransform(bcg, assay = "spatial", verbose = FALSE)

ntm_bcg <- SCTransform(ntm_bcg, assay = "spatial", verbose = FALSE)

dimensionality reduction and clustering

# saline
saline_analyzed <- RunPCA(saline, assay = "SCT", verbose = FALSE)
saline_analyzed <- FindNeighbors(saline_analyzed, reduction = "pca", dims = 1:5)
## Computing nearest neighbor graph
## Computing SNN
saline_analyzed <- FindClusters(saline_analyzed, verbose = FALSE, resolution = 0.8)
saline_analyzed <- RunUMAP(saline_analyzed, reduction = "pca", dims = 1:5)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 19:20:53 UMAP embedding parameters a = 0.9922 b = 1.112
## 19:20:53 Read 3049 rows and found 5 numeric columns
## 19:20:53 Using Annoy for neighbor search, n_neighbors = 30
## 19:20:53 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 19:20:53 Writing NN index file to temp file /var/folders/bx/5sy1zly51q9fkk63tcjztkr40000gq/T//RtmpiCls1V/file2acd94e3a12
## 19:20:53 Searching Annoy index using 1 thread, search_k = 3000
## 19:20:54 Annoy recall = 100%
## 19:20:55 Commencing smooth kNN distance calibration using 1 thread
## 19:20:56 Initializing from normalized Laplacian + noise
## 19:20:56 Commencing optimization for 500 epochs, with 116086 positive edges
## 19:21:01 Optimization finished
# ntm
ntm_analyzed <- RunPCA(ntm, assay = "SCT", verbose = FALSE)
ntm_analyzed <- FindNeighbors(ntm_analyzed, reduction = "pca", dims = 1:5)
## Computing nearest neighbor graph
## Computing SNN
ntm_analyzed <- FindClusters(ntm_analyzed, verbose = FALSE, resolution = 0.8)
ntm_analyzed <- RunUMAP(ntm_analyzed, reduction = "pca", dims = 1:5)
## 19:21:08 UMAP embedding parameters a = 0.9922 b = 1.112
## 19:21:08 Read 2174 rows and found 5 numeric columns
## 19:21:08 Using Annoy for neighbor search, n_neighbors = 30
## 19:21:08 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 19:21:08 Writing NN index file to temp file /var/folders/bx/5sy1zly51q9fkk63tcjztkr40000gq/T//RtmpiCls1V/file2acd6ec67883
## 19:21:08 Searching Annoy index using 1 thread, search_k = 3000
## 19:21:09 Annoy recall = 100%
## 19:21:09 Commencing smooth kNN distance calibration using 1 thread
## 19:21:10 Initializing from normalized Laplacian + noise
## 19:21:10 Commencing optimization for 500 epochs, with 84330 positive edges
## 19:21:15 Optimization finished
# bcg
bcg_analyzed <- RunPCA(bcg, assay = "SCT", verbose = FALSE)
bcg_analyzed <- FindNeighbors(bcg_analyzed, reduction = "pca", dims = 1:5)
## Computing nearest neighbor graph
## Computing SNN
bcg_analyzed <- FindClusters(bcg_analyzed, verbose = FALSE, resolution = 0.8)
bcg_analyzed <- RunUMAP(bcg_analyzed, reduction = "pca", dims = 1:5)
## 19:21:27 UMAP embedding parameters a = 0.9922 b = 1.112
## 19:21:27 Read 3318 rows and found 5 numeric columns
## 19:21:27 Using Annoy for neighbor search, n_neighbors = 30
## 19:21:27 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 19:21:27 Writing NN index file to temp file /var/folders/bx/5sy1zly51q9fkk63tcjztkr40000gq/T//RtmpiCls1V/file2acd28121946
## 19:21:27 Searching Annoy index using 1 thread, search_k = 3000
## 19:21:28 Annoy recall = 100%
## 19:21:29 Commencing smooth kNN distance calibration using 1 thread
## 19:21:30 Initializing from normalized Laplacian + noise
## 19:21:30 Commencing optimization for 500 epochs, with 126950 positive edges
## 19:21:35 Optimization finished
# ntm+bcg
ntm_bcg_analyzed <- RunPCA(ntm_bcg, assay = "SCT", verbose = FALSE)
ntm_bcg_analyzed <- FindNeighbors(ntm_bcg_analyzed, reduction = "pca", dims = 1:5)
## Computing nearest neighbor graph
## Computing SNN
ntm_bcg_analyzed <- FindClusters(ntm_bcg_analyzed, verbose = FALSE, resolution = 0.8)
ntm_bcg_analyzed <- RunUMAP(ntm_bcg_analyzed, reduction = "pca", dims = 1:5)
## 19:21:40 UMAP embedding parameters a = 0.9922 b = 1.112
## 19:21:40 Read 1947 rows and found 5 numeric columns
## 19:21:40 Using Annoy for neighbor search, n_neighbors = 30
## 19:21:40 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 19:21:40 Writing NN index file to temp file /var/folders/bx/5sy1zly51q9fkk63tcjztkr40000gq/T//RtmpiCls1V/file2acd3c395737
## 19:21:40 Searching Annoy index using 1 thread, search_k = 3000
## 19:21:41 Annoy recall = 100%
## 19:21:42 Commencing smooth kNN distance calibration using 1 thread
## 19:21:43 Initializing from normalized Laplacian + noise
## 19:21:43 Commencing optimization for 500 epochs, with 74996 positive edges
## 19:21:47 Optimization finished

labelling clusters using SingleR pipeline

# saline
saline_analyzed.sce <- as.SingleCellExperiment(saline_analyzed)

# ntm
ntm_analyzed.sce <- as.SingleCellExperiment(ntm_analyzed)

# bcg
bcg_analyzed.sce <- as.SingleCellExperiment(bcg_analyzed)

# ntm
ntm_bcg_analyzed.sce <- as.SingleCellExperiment(ntm_bcg_analyzed)

#create reference data

ref.data <- celldex::ImmGenData()
## snapshotDate(): 2021-10-19
## see ?celldex and browseVignettes('celldex') for documentation
## loading from cache
## see ?celldex and browseVignettes('celldex') for documentation
## loading from cache
# run singleR pipeline to find the cell types

## saline
cell_types_saline <- SingleR(test=saline_analyzed.sce, assay.type.test=1, 
    ref=ref.data, labels=ref.data$label.main)

## ntm
cell_types_ntm <- SingleR(test=ntm_analyzed.sce, assay.type.test=1, 
    ref=ref.data, labels=ref.data$label.main)

#bcg
cell_types_bcg <- SingleR(test=bcg_analyzed.sce, assay.type.test=1, 
    ref=ref.data, labels=ref.data$label.main)

#ntm+bcg
cell_types_ntm_bcg <- SingleR(test=ntm_bcg_analyzed.sce, assay.type.test=1, 
    ref=ref.data, labels=ref.data$label.main)

View cell types

# saline
table(cell_types_saline$labels)
## 
##           B cells                DC Endothelial cells  Epithelial cells 
##                58                41              1095               262 
##       Fibroblasts               ILC       Macrophages         Monocytes 
##               384                 1              1007                19 
##       Neutrophils     Stromal cells 
##                 2               180
# saline
table(cell_types_ntm$labels)
## 
##           B cells                DC Endothelial cells  Epithelial cells 
##               124                19               970                28 
##       Fibroblasts       Macrophages         Monocytes       Neutrophils 
##               234               477                 6                 4 
##     Stromal cells 
##               312
# saline
table(cell_types_bcg$labels)
## 
##           B cells                DC Endothelial cells  Epithelial cells 
##                48                43              1244               264 
##       Fibroblasts       Macrophages         Monocytes       Neutrophils 
##               437              1052                20                 2 
##     Stromal cells 
##               208
# saline
table(cell_types_ntm_bcg$labels)
## 
##           B cells                DC Endothelial cells  Epithelial cells 
##               192                93               578               100 
##       Fibroblasts       Macrophages         Monocytes     Stromal cells 
##               131               762                16                75

merging metadata of cell types to seurat object

# saline
saline_analyzed[["ref.data"]] <- cell_types_saline$labels

# saline
ntm_analyzed[["ref.data"]] <- cell_types_ntm$labels

# saline
bcg_analyzed[["ref.data"]] <- cell_types_bcg$labels

# saline
ntm_bcg_analyzed[["ref.data"]] <- cell_types_ntm_bcg$labels

my colors

my_cols <- c('B cells'='#F8766D','DC'='#39568CFF','Endothelial cells'='#CD9600','Epithelial cells'='#00C19A','Fibroblasts'= "#C77CFF",
   "ILC" = "grey", 'Macrophages'='#00A9FF','Monocytes'='#0CB702', "Neutrophils" = "Darkred", 'Stromal cells'='#FF61CC')

Create cell chat object and merge

saline

saline_cellchat <- createCellChat(object = saline_analyzed, group.by = "ref.data")
## Create a CellChat object from a Seurat object
## The `data` slot in the default assay is used. The default assay is SCT 
## The `meta.data` slot in the Seurat object is used as cell meta information
## Set cell identities for the new CellChat object
## The cell groups used for CellChat analysis are  B cells DC Endothelial cells Epithelial cells Fibroblasts ILC Macrophages Monocytes Neutrophils Stromal cells
CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)

CellChatDB.use <- CellChatDB
saline_cellchat@DB <- CellChatDB.use


saline_cellchat <- subsetData(saline_cellchat)
saline_cellchat <- identifyOverExpressedGenes(saline_cellchat)
saline_cellchat <- identifyOverExpressedInteractions(saline_cellchat)


saline_cellchat <- computeCommunProb(saline_cellchat)
saline_cellchat <- filterCommunication(saline_cellchat, min.cells = 10)
## The cell-cell communication related with the following cell groups are excluded due to the few number of cells:  ILC Neutrophils
saline_cellchat <- computeCommunProbPathway(saline_cellchat)
saline_cellchat <- aggregateNet(saline_cellchat)

groupSize <- as.numeric(table(saline_cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(saline_cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions", color.use = my_cols)
netVisual_circle(saline_cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength", color.use = my_cols)

# See interaction of each cell type with the another

mat1 <- saline_cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat1)) {
  mat2 <- matrix(0, nrow = nrow(mat1), ncol = ncol(mat1), dimnames = dimnames(mat1))
  mat2[i, ] <- mat1[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat1), title.name = rownames(mat1)[i], color.use = my_cols)}

Visualize signalling pathways in each cell

# Identify signaling pathways showing significant communications 
saline_cellchat@netP$pathways
##  [1] "COLLAGEN"    "APP"         "MIF"         "THBS"        "VEGF"       
##  [6] "FN1"         "LAMININ"     "MHC-II"      "ICAM"        "ITGAL-ITGB2"
## [11] "COMPLEMENT"  "PECAM1"      "CD52"        "CD22"        "CD45"       
## [16] "JAM"         "PARs"        "CHEMERIN"    "SEMA4"       "CXCL"       
## [21] "EPHA"        "PDGF"        "GAS"         "HSPG"        "FGF"        
## [26] "SEMA3"       "AGRN"        "PTPRM"       "CDH5"        "SEMA7"      
## [31] "NOTCH"       "SEMA6"       "IL16"
# Evaluate one of the pathways
pathways.show <- c("MHC-II") 
vertex.receiver = seq(2,5) # a numeric vector. 
netVisual_aggregate(saline_cellchat, signaling = pathways.show,  vertex.receiver = vertex.receiver, color.use = my_cols)

# Chord diagram
par(mfrow=c(1,1))
netVisual_aggregate(saline_cellchat, signaling = pathways.show, layout = "chord", color.use = my_cols)

# Heatmap
par(mfrow=c(1,1))
netVisual_heatmap(saline_cellchat, signaling = pathways.show, color.heatmap = "Reds", color.use = my_cols)
## Do heatmap based on a single object

Compute the contribution of each ligand-receptor pair to the overall signaling pathway

netAnalysis_contribution(saline_cellchat, signaling = pathways.show)

pairLR.MHC <- extractEnrichedLR(saline_cellchat, signaling = pathways.show, geneLR.return = FALSE)
LR.show <- pairLR.MHC[1,] # show one ligand-receptor pair

# Hierarchy plot
vertex.receiver = seq(1,4) # a numeric vector
netVisual_individual(saline_cellchat, signaling = pathways.show,  pairLR.use = LR.show, vertex.receiver = vertex.receiver, color.use = my_cols)

## [[1]]

Visualize cell-cell communication mediated by multiple signaling pathways

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(saline_cellchat, sources.use = 1, targets.use = c(2:8), remove.isolate = FALSE, font.size = 14)
## Comparing communications on a single object

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB
netVisual_chord_gene(saline_cellchat, sources.use = 1, targets.use = c(2:8), lab.cex = 0.7,legend.pos.y = 30, color.use = my_cols)

netVisual_chord_gene(saline_cellchat, sources.use = 7, targets.use = c(1:8), lab.cex = 0.7,legend.pos.y = 30, color.use = my_cols)

Compute the network centrality scores

saline_cellchat <- netAnalysis_computeCentrality(saline_cellchat, slot.name = "netP") 

# the slot 'netP' means the inferred intercellular communication network of signaling pathways

# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
netAnalysis_signalingRole_network(saline_cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10, color.use = my_cols)

Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways

ht1 <- netAnalysis_signalingRole_heatmap(saline_cellchat, pattern = "outgoing", color.use = my_cols)
ht2 <- netAnalysis_signalingRole_heatmap(saline_cellchat, pattern = "incoming", color.use = my_cols)
ht1 + ht2

Identify and visualize communication patterns (outgoing signal)

selectK(saline_cellchat, pattern = "outgoing")

nPatterns = 3
saline_cellchat <- identifyCommunicationPatterns(saline_cellchat, pattern = "outgoing", k = nPatterns)

# river plot
netAnalysis_river(saline_cellchat, pattern = "outgoing", color.use = my_cols)
## Please make sure you have load `library(ggalluvial)` when running this function

# dot plot
netAnalysis_dot(saline_cellchat, pattern = "outgoing", color.use = my_cols)

Identify and visualize communication patterns (incoming signal)

selectK(saline_cellchat, pattern = "incoming")

nPatterns = 3
saline_cellchat <- identifyCommunicationPatterns(saline_cellchat, pattern = "incoming", k = nPatterns)

# river plot
netAnalysis_river(saline_cellchat, pattern = "incoming", color.use = my_cols)
## Please make sure you have load `library(ggalluvial)` when running this function

# dot plot
netAnalysis_dot(saline_cellchat, pattern = "incoming", color.use = my_cols)

saline_cellchat <- computeNetSimilarity(saline_cellchat, type = "functional")
saline_cellchat <- netEmbedding(saline_cellchat, type = "functional")
## Manifold learning of the signaling networks for a single dataset
#> Manifold learning of the signaling networks for a single dataset
saline_cellchat <- netClustering(saline_cellchat, type = "functional")
## Classification learning of the signaling networks for a single dataset
#> Classification learning of the signaling networks for a single dataset
# Visualization in 2D-space
netVisual_embedding(saline_cellchat, type = "functional", label.size = 3.5, color.use = my_cols)

NTM

my colors

my_cols <- c('B cells'='#F8766D','DC'='#39568CFF','Endothelial cells'='#CD9600','Epithelial cells'='#00C19A','Fibroblasts'= "#C77CFF",
    'Macrophages'='#00A9FF','Monocytes'='#0CB702', "Neutrophils" = "Darkred", 'Stromal cells'='#FF61CC')
ntm_cellchat <- createCellChat(object = ntm_analyzed, group.by = "ref.data")
## Create a CellChat object from a Seurat object
## The `data` slot in the default assay is used. The default assay is SCT 
## The `meta.data` slot in the Seurat object is used as cell meta information
## Set cell identities for the new CellChat object
## The cell groups used for CellChat analysis are  B cells DC Endothelial cells Epithelial cells Fibroblasts Macrophages Monocytes Neutrophils Stromal cells
CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)

CellChatDB.use <- CellChatDB
ntm_cellchat@DB <- CellChatDB.use


ntm_cellchat <- subsetData(ntm_cellchat)
ntm_cellchat <- identifyOverExpressedGenes(ntm_cellchat)
ntm_cellchat <- identifyOverExpressedInteractions(ntm_cellchat)


ntm_cellchat <- computeCommunProb(ntm_cellchat)
ntm_cellchat <- filterCommunication(ntm_cellchat, min.cells = 10)
## The cell-cell communication related with the following cell groups are excluded due to the few number of cells:  Monocytes Neutrophils
ntm_cellchat <- computeCommunProbPathway(ntm_cellchat)
ntm_cellchat <- aggregateNet(ntm_cellchat)

groupSize <- as.numeric(table(ntm_cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(ntm_cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions", color.use = my_cols)
netVisual_circle(ntm_cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength", color.use = my_cols)

# See interaction of each cell type with the another

mat1 <- ntm_cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat1)) {
  mat2 <- matrix(0, nrow = nrow(mat1), ncol = ncol(mat1), dimnames = dimnames(mat1))
  mat2[i, ] <- mat1[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat1), title.name = rownames(mat1)[i], color.use = my_cols)
}

Visualize signalling pathways in each cell

# Identify signaling pathways showing significant communications 
ntm_cellchat@netP$pathways
##  [1] "COLLAGEN"    "APP"         "THBS"        "FN1"         "MIF"        
##  [6] "COMPLEMENT"  "GALECTIN"    "MHC-II"      "VEGF"        "LAMININ"    
## [11] "SEMA4"       "ICAM"        "CXCL"        "GAS"         "CCL"        
## [16] "ITGAL-ITGB2" "JAM"         "PECAM1"      "CD22"        "CD45"       
## [21] "PDGF"        "CD52"        "TENASCIN"    "THY1"        "TGFb"       
## [26] "FGF"         "SEMA3"       "CHEMERIN"    "PD-L1"       "LT"         
## [31] "SEMA7"       "SPP1"        "NOTCH"       "GRN"         "SEMA6"      
## [36] "ESAM"        "IL16"        "EPHB"        "PARs"        "CDH"        
## [41] "ncWNT"       "IGF"         "CD39"        "APRIL"       "ICOS"       
## [46] "PDL2"        "NRG"         "XCR"         "IL1"
# Evaluate one of the pathways
pathways.show <- c("MHC-II") 
vertex.receiver = seq(2,5) # a numeric vector. 
netVisual_aggregate(ntm_cellchat, signaling = pathways.show,  vertex.receiver = vertex.receiver, color.use = my_cols)

# Chord diagram
par(mfrow=c(1,1))
netVisual_aggregate(ntm_cellchat, signaling = pathways.show, layout = "chord", color.use = my_cols)

# Heatmap
par(mfrow=c(1,1))
netVisual_heatmap(ntm_cellchat, signaling = pathways.show, color.heatmap = "Reds", color.use = my_cols)
## Do heatmap based on a single object

Compute the contribution of each ligand-receptor pair to the overall signaling pathway

netAnalysis_contribution(ntm_cellchat, signaling = pathways.show)

pairLR.MHC <- extractEnrichedLR(ntm_cellchat, signaling = pathways.show, geneLR.return = FALSE)
LR.show <- pairLR.MHC[1,] # show one ligand-receptor pair

# Hierarchy plot
vertex.receiver = seq(1,4) # a numeric vector
netVisual_individual(ntm_cellchat, signaling = pathways.show,  pairLR.use = LR.show, vertex.receiver = vertex.receiver, color.use = my_cols)

## [[1]]

Visualize cell-cell communication mediated by multiple signaling pathways

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(ntm_cellchat, sources.use = 1, targets.use = c(2:8), remove.isolate = FALSE, font.size = 14)
## Comparing communications on a single object

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB
netVisual_chord_gene(ntm_cellchat, sources.use = 1, targets.use = c(2:8), lab.cex = 0.7,legend.pos.y = 30, color.use = my_cols)

Compute the network centrality scores

ntm_cellchat <- netAnalysis_computeCentrality(ntm_cellchat, slot.name = "netP") 

# the slot 'netP' means the inferred intercellular communication network of signaling pathways

# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
netAnalysis_signalingRole_network(ntm_cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10, color.use = my_cols)

Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways

ht1 <- netAnalysis_signalingRole_heatmap(ntm_cellchat, pattern = "outgoing", color.use = my_cols)
ht2 <- netAnalysis_signalingRole_heatmap(ntm_cellchat, pattern = "incoming", color.use = my_cols)
ht1 + ht2

Identify and visualize communication patterns (outgoing signal)

selectK(ntm_cellchat, pattern = "outgoing")

nPatterns = 3
ntm_cellchat <- identifyCommunicationPatterns(ntm_cellchat, pattern = "outgoing", k = nPatterns)

# river plot
netAnalysis_river(ntm_cellchat, pattern = "outgoing", color.use = my_cols)
## Please make sure you have load `library(ggalluvial)` when running this function

# dot plot
netAnalysis_dot(ntm_cellchat, pattern = "outgoing", color.use = my_cols)

Identify and visualize communication patterns (incoming signal)

selectK(ntm_cellchat, pattern = "incoming")

nPatterns = 3
ntm_cellchat <- identifyCommunicationPatterns(ntm_cellchat, pattern = "incoming", k = nPatterns)

# river plot
netAnalysis_river(ntm_cellchat, pattern = "incoming", color.use = my_cols)
## Please make sure you have load `library(ggalluvial)` when running this function

# dot plot
netAnalysis_dot(ntm_cellchat, pattern = "incoming", color.use = my_cols)

ntm_cellchat <- computeNetSimilarity(ntm_cellchat, type = "functional")
ntm_cellchat <- netEmbedding(ntm_cellchat, type = "functional")
## Manifold learning of the signaling networks for a single dataset
#> Manifold learning of the signaling networks for a single dataset
ntm_cellchat <- netClustering(ntm_cellchat, type = "functional")
## Classification learning of the signaling networks for a single dataset
#> Classification learning of the signaling networks for a single dataset
# Visualization in 2D-space
netVisual_embedding(ntm_cellchat, type = "functional", label.size = 3.5, color.use = my_cols)

BCG

bcg_cellchat <- createCellChat(object = bcg_analyzed, group.by = "ref.data")
## Create a CellChat object from a Seurat object
## The `data` slot in the default assay is used. The default assay is SCT 
## The `meta.data` slot in the Seurat object is used as cell meta information
## Set cell identities for the new CellChat object
## The cell groups used for CellChat analysis are  B cells DC Endothelial cells Epithelial cells Fibroblasts Macrophages Monocytes Neutrophils Stromal cells
CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)

CellChatDB.use <- CellChatDB
bcg_cellchat@DB <- CellChatDB.use


bcg_cellchat <- subsetData(bcg_cellchat)
bcg_cellchat <- identifyOverExpressedGenes(bcg_cellchat)
bcg_cellchat <- identifyOverExpressedInteractions(bcg_cellchat)


bcg_cellchat <- computeCommunProb(bcg_cellchat)
bcg_cellchat <- filterCommunication(bcg_cellchat, min.cells = 10)
## The cell-cell communication related with the following cell groups are excluded due to the few number of cells:  Neutrophils
bcg_cellchat <- computeCommunProbPathway(bcg_cellchat)
bcg_cellchat <- aggregateNet(bcg_cellchat)

groupSize <- as.numeric(table(bcg_cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(bcg_cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions", color.use = my_cols)
netVisual_circle(bcg_cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength", color.use = my_cols)

# See interaction of each cell type with the another

mat1 <- bcg_cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat1)) {
  mat2 <- matrix(0, nrow = nrow(mat1), ncol = ncol(mat1), dimnames = dimnames(mat1))
  mat2[i, ] <- mat1[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat1), title.name = rownames(mat1)[i], color.use = my_cols)
}

Visualize signalling pathways in each cell

# Identify signaling pathways showing significant communications 
bcg_cellchat@netP$pathways
##  [1] "APP"         "COLLAGEN"    "GALECTIN"    "THBS"        "MIF"        
##  [6] "LAMININ"     "VEGF"        "MHC-II"      "FN1"         "ICAM"       
## [11] "COMPLEMENT"  "JAM"         "ITGAL-ITGB2" "SPP1"        "PECAM1"     
## [16] "CHEMERIN"    "CD22"        "CD45"        "CD52"        "GAS"        
## [21] "SEMA3"       "CXCL"        "SEMA4"       "FGF"         "SEMA7"      
## [26] "IL16"        "PDGF"
# Evaluate one of the pathways
pathways.show <- c("MHC-II") 
vertex.receiver = seq(2,5) # a numeric vector. 
netVisual_aggregate(bcg_cellchat, signaling = pathways.show,  vertex.receiver = vertex.receiver, color.use = my_cols)

# Chord diagram
par(mfrow=c(1,1))
netVisual_aggregate(bcg_cellchat, signaling = pathways.show, layout = "chord", color.use = my_cols)

# Heatmap
par(mfrow=c(1,1))
netVisual_heatmap(bcg_cellchat, signaling = pathways.show, color.heatmap = "Reds", color.use = my_cols)
## Do heatmap based on a single object

Compute the contribution of each ligand-receptor pair to the overall signaling pathway

netAnalysis_contribution(bcg_cellchat, signaling = pathways.show)

pairLR.MHC <- extractEnrichedLR(bcg_cellchat, signaling = pathways.show, geneLR.return = FALSE)
LR.show <- pairLR.MHC[1,] # show one ligand-receptor pair

# Hierarchy plot
vertex.receiver = seq(1,4) # a numeric vector
netVisual_individual(bcg_cellchat, signaling = pathways.show,  pairLR.use = LR.show, vertex.receiver = vertex.receiver, color.use = my_cols)

## [[1]]

Visualize cell-cell communication mediated by multiple signaling pathways

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(bcg_cellchat, sources.use = 1, targets.use = c(2:8), remove.isolate = FALSE, font.size = 14)
## Comparing communications on a single object

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB
netVisual_chord_gene(bcg_cellchat, sources.use = 1, targets.use = c(2:8), lab.cex = 0.7,legend.pos.y = 30, color.use = my_cols)

Compute the network centrality scores

bcg_cellchat <- netAnalysis_computeCentrality(bcg_cellchat, slot.name = "netP") 

# the slot 'netP' means the inferred intercellular communication network of signaling pathways

# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
netAnalysis_signalingRole_network(bcg_cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10, color.use = my_cols)

Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways

ht1 <- netAnalysis_signalingRole_heatmap(bcg_cellchat, pattern = "outgoing", color.use = my_cols)
ht2 <- netAnalysis_signalingRole_heatmap(bcg_cellchat, pattern = "incoming", color.use = my_cols)
ht1 + ht2

Identify and visualize communication patterns (outgoing signal)

selectK(bcg_cellchat, pattern = "outgoing")

nPatterns = 3
bcg_cellchat <- identifyCommunicationPatterns(bcg_cellchat, pattern = "outgoing", k = nPatterns)

# river plot
netAnalysis_river(bcg_cellchat, pattern = "outgoing", color.use = my_cols)
## Please make sure you have load `library(ggalluvial)` when running this function

# dot plot
netAnalysis_dot(bcg_cellchat, pattern = "outgoing", color.use = my_cols)

Identify and visualize communication patterns (incoming signal)

selectK(bcg_cellchat, pattern = "incoming")

nPatterns = 3
bcg_cellchat <- identifyCommunicationPatterns(bcg_cellchat, pattern = "incoming", k = nPatterns)

# river plot
netAnalysis_river(bcg_cellchat, pattern = "incoming", color.use = my_cols)
## Please make sure you have load `library(ggalluvial)` when running this function

# dot plot
netAnalysis_dot(bcg_cellchat, pattern = "incoming", color.use = my_cols)

bcg_cellchat <- computeNetSimilarity(bcg_cellchat, type = "functional")
bcg_cellchat <- netEmbedding(bcg_cellchat, type = "functional")
## Manifold learning of the signaling networks for a single dataset
#> Manifold learning of the signaling networks for a single dataset
bcg_cellchat <- netClustering(bcg_cellchat, type = "functional")
## Classification learning of the signaling networks for a single dataset
#> Classification learning of the signaling networks for a single dataset
# Visualization in 2D-space
netVisual_embedding(bcg_cellchat, type = "functional", label.size = 3.5, color.use = my_cols)

NTM + BCG

my colors

my_cols <- c('B cells'='#F8766D','DC'='#39568CFF','Endothelial cells'='#CD9600','Epithelial cells'='#00C19A','Fibroblasts'= "#C77CFF",
    'Macrophages'='#00A9FF','Monocytes'='#0CB702', 'Stromal cells'='#FF61CC')
ntm_bcg_cellchat <- createCellChat(object = ntm_bcg_analyzed, group.by = "ref.data")
## Create a CellChat object from a Seurat object
## The `data` slot in the default assay is used. The default assay is SCT 
## The `meta.data` slot in the Seurat object is used as cell meta information
## Set cell identities for the new CellChat object
## The cell groups used for CellChat analysis are  B cells DC Endothelial cells Epithelial cells Fibroblasts Macrophages Monocytes Stromal cells
CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)

CellChatDB.use <- CellChatDB
ntm_bcg_cellchat@DB <- CellChatDB.use


ntm_bcg_cellchat <- subsetData(ntm_bcg_cellchat)
ntm_bcg_cellchat <- identifyOverExpressedGenes(ntm_bcg_cellchat)
ntm_bcg_cellchat <- identifyOverExpressedInteractions(ntm_bcg_cellchat)


ntm_bcg_cellchat <- computeCommunProb(ntm_bcg_cellchat)
ntm_bcg_cellchat <- filterCommunication(ntm_bcg_cellchat, min.cells = 10)
ntm_bcg_cellchat <- computeCommunProbPathway(ntm_bcg_cellchat)
ntm_bcg_cellchat <- aggregateNet(ntm_bcg_cellchat)

groupSize <- as.numeric(table(ntm_bcg_cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(ntm_bcg_cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions", color.use = my_cols)
netVisual_circle(ntm_bcg_cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength", color.use = my_cols)

# See interaction of each cell type with the another

mat1 <- ntm_bcg_cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat1)) {
  mat2 <- matrix(0, nrow = nrow(mat1), ncol = ncol(mat1), dimnames = dimnames(mat1))
  mat2[i, ] <- mat1[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat1), title.name = rownames(mat1)[i], color.use = my_cols)
}

Visualize signalling pathways in each cell

# Identify signaling pathways showing significant communications 
ntm_bcg_cellchat@netP$pathways
##  [1] "COLLAGEN"    "APP"         "THBS"        "MIF"         "GALECTIN"   
##  [6] "MHC-II"      "COMPLEMENT"  "LAMININ"     "GAS"         "FN1"        
## [11] "JAM"         "SEMA4"       "VEGF"        "SEMA3"       "CXCL"       
## [16] "ICAM"        "CCL"         "SPP1"        "ITGAL-ITGB2" "PECAM1"     
## [21] "CD22"        "CD45"        "CD52"        "LT"          "THY1"       
## [26] "CHEMERIN"    "TENASCIN"    "FGF"         "SEMA7"       "TNF"        
## [31] "WNT"         "PD-L1"       "PARs"        "PDGF"        "CADM"       
## [36] "GRN"         "IL16"        "PROS"        "ESAM"        "TWEAK"      
## [41] "CDH"         "ALCAM"       "CD6"         "TGFb"        "ncWNT"      
## [46] "ICOS"        "CSF"         "AGRN"        "APRIL"       "PDL2"       
## [51] "CDH5"
# Evaluate one of the pathways
pathways.show <- c("MHC-II") 
vertex.receiver = seq(2,5) # a numeric vector. 
netVisual_aggregate(ntm_bcg_cellchat, signaling = pathways.show,  vertex.receiver = vertex.receiver, color.use = my_cols)

# Chord diagram
par(mfrow=c(1,1))
netVisual_aggregate(ntm_bcg_cellchat, signaling = pathways.show, layout = "chord", color.use = my_cols)

# Heatmap
par(mfrow=c(1,1))
netVisual_heatmap(ntm_bcg_cellchat, signaling = pathways.show, color.heatmap = "Reds", color.use = my_cols)
## Do heatmap based on a single object

Compute the contribution of each ligand-receptor pair to the overall signaling pathway

netAnalysis_contribution(ntm_bcg_cellchat, signaling = pathways.show)

pairLR.MHC <- extractEnrichedLR(ntm_bcg_cellchat, signaling = pathways.show, geneLR.return = FALSE)
LR.show <- pairLR.MHC[1,] # show one ligand-receptor pair

# Hierarchy plot
vertex.receiver = seq(1,4) # a numeric vector
netVisual_individual(ntm_bcg_cellchat, signaling = pathways.show,  pairLR.use = LR.show, vertex.receiver = vertex.receiver, color.use = my_cols)

## [[1]]

Visualize cell-cell communication mediated by multiple signaling pathways

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(ntm_bcg_cellchat, sources.use = 1, targets.use = c(2:8), remove.isolate = FALSE, font.size = 14)
## Comparing communications on a single object

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB
netVisual_chord_gene(ntm_bcg_cellchat, sources.use = 1, targets.use = c(2:8), lab.cex = 0.7, legend.pos.y = 30, color.use = my_cols)

Compute the network centrality scores

ntm_bcg_cellchat <- netAnalysis_computeCentrality(ntm_bcg_cellchat, slot.name = "netP") 

# the slot 'netP' means the inferred intercellular communication network of signaling pathways

# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
netAnalysis_signalingRole_network(ntm_bcg_cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10, color.use = my_cols)

Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways

ht1 <- netAnalysis_signalingRole_heatmap(ntm_bcg_cellchat, pattern = "outgoing", color.use = my_cols)
ht2 <- netAnalysis_signalingRole_heatmap(ntm_bcg_cellchat, pattern = "incoming", color.use = my_cols)
ht1 + ht2

Identify and visualize communication patterns (outgoing signal)

selectK(ntm_bcg_cellchat, pattern = "outgoing")

nPatterns = 3
ntm_bcg_cellchat <- identifyCommunicationPatterns(ntm_bcg_cellchat, pattern = "outgoing", k = nPatterns)

# river plot
netAnalysis_river(ntm_bcg_cellchat, pattern = "outgoing", color.use = my_cols)
## Please make sure you have load `library(ggalluvial)` when running this function

# dot plot
netAnalysis_dot(ntm_bcg_cellchat, pattern = "outgoing", color.use = my_cols)

Identify and visualize communication patterns (incoming signal)

selectK(ntm_bcg_cellchat, pattern = "incoming")

nPatterns = 3
ntm_bcg_cellchat <- identifyCommunicationPatterns(ntm_bcg_cellchat, pattern = "incoming", k = nPatterns)

# river plot
netAnalysis_river(ntm_bcg_cellchat, pattern = "incoming", color.use = my_cols)
## Please make sure you have load `library(ggalluvial)` when running this function

# dot plot
netAnalysis_dot(ntm_bcg_cellchat, pattern = "incoming", color.use = my_cols)

ntm_bcg_cellchat <- computeNetSimilarity(ntm_bcg_cellchat, type = "functional")
ntm_bcg_cellchat <- netEmbedding(ntm_bcg_cellchat, type = "functional")
## Manifold learning of the signaling networks for a single dataset
#> Manifold learning of the signaling networks for a single dataset
ntm_bcg_cellchat <- netClustering(ntm_bcg_cellchat, type = "functional")
## Classification learning of the signaling networks for a single dataset
#> Classification learning of the signaling networks for a single dataset
# Visualization in 2D-space
netVisual_embedding(ntm_bcg_cellchat, type = "functional", label.size = 3.5, color.use = my_cols)

merge objects

object.list <- list(saline_group = saline_cellchat, ntm_group = ntm_cellchat, bcg_group = bcg_cellchat, ntm_bcg_group = ntm_bcg_cellchat)
cellchat <- mergeCellChat(object.list, add.names = names(object.list), cell.prefix = TRUE)
## The cell barcodes in merged 'meta' is  AAACAAGTATCTCCCA-1 AAACAGAGCGACTCCT-1 AAACATTTCCCGGATT-1 AAACCCGAACGAAATC-1 AAACCGGGTAGGTACC-1 AAACCGTTCGTCCAGG-1
## Merge the following slots: 'data.signaling','net', 'netP','meta', 'idents', 'var.features' , 'DB', and 'LR'.

compare interactions

my_cols_group <- c("red2", "green4", "grey59", "dodgerblue3")

gg1 <- compareInteractions(cellchat, show.legend = F, group = c(1,2,3,4), color.use = my_cols_group)
gg2 <- compareInteractions(cellchat, show.legend = F, group = c(1,2,3,4), measure = "weight", color.use = my_cols_group)
gg1 + gg2